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ABSTRACT 


This study involves the development of numerical modelling in dilute and dense 
spray combustion. We discuss several issues concerning the computational effi- 
ciency in the stochastic particle tracking method as well as the improvement of 
physical submodels of turbulence, combustion, vaporization, swirling effects, and 
dense spray effects. The governing gas-phase equations in Eulerian coordinate are 
solved by a time-marching multiple pressure correction procedure based on the 
operator-splitting technique. The droplet-phase equations in Lagrangian coordi- 
nate are solved by a stochastic discrete droplet technique. The k — e model is used 
to characterize the time and length scales of the gas phase turbulence for droplet 
dispersions and droplet /turbulence interactions. To improve the computational ef- 
ficiency in the stochastic tracking calculations, we implement a dispersion width 
transport model which can account for turbulent dispersion within each computa- 
tional parcel. The testings of this model confirm the capability of accuratly repre- 
senting dispersion in nearly-homogeneous and inhomogeneous turbulent flows with 
im proved efficiency over the delta function stochastic separated flow(SSF) model. 
To account for the dense spray effects, we have employed an existing drop collision 
and coalescence model and a Taylor analogy breakup(TAB) model. These models 
were incorporated into a state-of-the-art multiphase all-speed transient flow solution 
procedure. A sequence of validation cases involving non-evaporating, evaporating, 
burning, dilute and dense spray cases are included. The research tasks concerning 
the development of multidimensional group particle tracking method and particle 
wall-boundary condition are separately documented in Appendix E. 
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NOMENCLATURE 


B m : mass transfer number 

Bt : heat transfer number 

C p : heat capacity of air 

Cp t d : heat capacity of droplet 

Cq : drag coefficent 

dp : droplet diameter 

D : diffusion coefficient in Fick’s law 

E : mean specific internal energy 

Ep : specific internal energy of particle 

/ : mixture fraction 

Fpi : particle drag force 

g : gravity 

h : enthalpy of gas phase 

h / : fuel vapor enthalpy 

K : undersampling correction factor 

k : turbulent kinetic energy 

L : latent heat 

Nt : number of droplets for each 

computational particle p 
Np : n um ber of computational particles 

P : probability 

p : mean pressure 

Pr t : turbulent Prandtl number 

Re p : particle Reynold’s number 

r p : droplet radius 

Ui : instantaneous velocity for gases 

Vi : instantaneos velocity for droplets 

S : source terms 

Sc t : turbulent Schmidt number 

T : gas temperature 


IV 



t : 

time 

T d : 

droplet temperature 

x : 

coordinate in the streamwise direction 

Y : 

mass fraction 

y : 

coordinate normal to the streamwise 
direction 

Greek 

Symbol 

P : 

density of gases 

/>/ : 

fuel vapor density 

/>«* : 

droplet density 

A* : 

viscosity 

A*< : 

eddy viscosity 


instantaneous scalar for gases 

€ : 

turbulent energy dissipation rate 


instantaneous mixture fraction 

r : 

particle relaxation time 

a : 

standard deviation of pdf 

ct 2 : 

variance of pdf 

Subscripts 

/: 

fuel 

9 ■ 

gas phase 

i : 

time step index in an eddy 

k : 

eddy index 

l : 

liquid phase 

P ■ 

particle or parcel 

rms : 

root mean square 

t : 

turbulent 

Supercripts 

— : 

density-averaged 

f . 

fluctuating 
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1. INTRODUCTION 

1.1 General 

There have been a number of research efforts[l-6] towards the development of 
numerical and physical models for spray combustion. Many aspects of sprays includ- 
ing fuel properties of droplets [1], multicomponent nature of fuel [2], evaporation 
models [3] have been studied and excellent reviews on analysis and measurements 
of sprays have also been given in Refs. [4,5,6]. These studies are motivated by 
the need for better understanding of multi-phase turbulent combustion processes 
as well as the demand for improving performance, stability, and emission control in 
industrial furnaces and propulsive systems such as gas turbine, ramjet engines, and 
space shuttle main engines. 

The prediction of the local flow properties of spray flames requires the solu- 
tions of multi-phase dynamics, and accounts for complex interactions between the 
dispersed droplets and the continuous gas-phase flows. Various approaches have 
been suggested to model the interphase transport phenomena. The methodologies 
for the spray combustion computations are largely classified as the discrete droplet 
model, the statistical droplet model, and the two-fluid continuum model. Com- 
parative performances for three approaches are well summarized in Ref. [6]. Among 
three models, the discrete droplet model has gained wide acceptance due to its com- 
putaional efficiency, the flexibility in handling poly-disperse spray, the convenient 
interphase coupling, and the elimination of numerical diffusion. With Eulerian- 
Lagrangian formulations in multi-phase flows, the stochastic separated flow(SSF) 
approach [5] categorized in the discrete droplet model is usually employed to ac- 
count for the turbulence effects on interphase transport. In the present stochastic 
separated flow model, the mathematical formulation of the two-phase flow and com- 
bustion processes comprises the Eulerian conservation equation for the gas phase 
and the Lagrangian equations for the fuel droplets. The link between two phases 
is mathematically expressed in terms of liquid/gas-phase interaction source terms 
in the gas-phase equations. The governing gas-phase equations in Eulerian coor- 
dinate are solved by a time-marching multiple pressure correction procedure based 
on the operator-splitting technique. The droplet-phase equations in Lagrangian 
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coordinate are solved by a stochastic discrete droplet technique. The k c model 
is used to characterize the time and length scales of the gas phase turbulence for 
droplet dispersions and droplet/turbulence interactions. The present vaporization 
model includes the effects of variable thermophysical properties, non-umtary Lewis 
number in the gas-film, the Stefan flow effect, and the effect of internal circulation 

and transient liquid heating. 

This study is mainly motivated to improve the physical submodels as well as 
to enhance the prediction capabilitiy over a wider range of fuel spray conditions 
and combustor geometries. In the following subsections, we adress several issues 
concerning the computational efficiency in the stochastic particle tracking method 
as well as the improvement of physical submodels of turbulence, combustion, va- 
porization, swirling effects, and dense spray effects. 

1.2 Turbulent Particle Dispersion 

In the stochastic separated flow(SSF) approach, each computational parcel 
represents a collection of liquid droplets having the same droplet characteristics 
and a random sampling technique is entailed for instantaneous gas flow properties 
based on a specified turbulence model and the resulting fluctuations are used in the 
droplet-phase Lagrangian computations for the droplet tracking. This stochastic 
process requires a large number of computational particles to produce satisfactory 
dispersion distributions even for rather dilute sprays. To circumvent this deficiency, 
Litchford and Jeng(7] proposed the dispersion width(group) approach in which each 
computational parcel represents a group of physical particles with a probability 
density distribution. This dispersion width model can account for the turbulent 
droplet dispersion within each group. Each group( width) grows due to the turbu- 
lent dispersion of droplets when the computational parcel travels in the Lagrangian 
coordinates. The mean position of each group, determined from a deterministic or 
stochastic Lagrangian tracking, is taken to represent the mean of its correspond- 
ing probability density function(PDF). The variance of each PDF is represented 
by a statistical mean-squared dispersion which depends on prior eddy interactions. 
Potential advantages of this method is to reduce the number of computational par- 
ticles which represent the spray dynamics and to obtain grid-independent solutions 
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for two-phase flows regardless of grid-refinement in the underlying Eulerian gas-flow 
calculations. Other advantages may include better representations of "group” evap- 
oration or "group” combustion models. The dispersion group model to be presented 
in this report is basically similar to the approach of Litchford and Jeng. However, 
the present procedure is somewhat different from the method proposed by Litch- 
ford and Jeng[7], in which the calculation of dispersion- related parameters needs 
summation of the entire time history. Furthermore, the present procedure is easy 
to program and requires less computer memory. To evaluate the present disper- 
sion width transport model and to calibrate the stochastic simulation of particle- 
turbulence interactions, the computations were performed for the particle dispersion 
in nearly-homogeneous turbulence and a particle laden round jet in inhomogeneous 
turbulence. The present dispersion width transport model has successfully demon- 
strated the capability of accuratly representing dispersion in nearly-homogeneous 
and inhomogeneous turbulent flows with improved efficiency over the delta function 

SSF model. 

1.3 Dilute Spray Models 

In the dilute spray combustion models, the stochastic separated flow model is 
employed to account for the turbulent droplet dispersion, turbulence is represented 
by the k — e model, and the combustion processes involves an irreversible one-step 
reaction at an infinite rate. The turbulent fluctuations on the mixture properties are 
introduced by the probability density function(pdf) approach. The centrifugal force 
terms associated with the swirl effects are also included in the gas-phase/droplet- 
phase equations. In the study, we evaluate the solution procedure and the physi- 
cal submodels of turbulence, combustion, vaporization, swirling effects, and initial 
spray distributions. The present numerical model for the multi-phase turbulent 
reacting flows has been tested by applying it to predict the local flow properties in 
two axisymmetric, confined, swirling spray-combusting flows[36,37). Example prob- 
lems include the liquid-fuel combusting flows with a hollow cone spray and with a 
rotating cup atomizer. Special emphasis is given to the influence of the spray initial 
conditions and the inlet swirl strength which characterize the spray vaporization 
and the turbulent mixing. Two swirl numbers are considered to investigate the 
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influence of swirl on the droplet evaporation and trajectories, and the effects of 
droplet/ turbulence interactions in flow properties. The predictive capabilities of 
the present procedure have been demonstrated by comparisons with experimen- 
tal data. The present numerical procedure correctly predicts the general features 
of spray-combustion flows and yields the qualitative agreement with experimental 
data. However, quantitative differences exist especially at near-burner locations, at 
near-wall regions, and along the combustion chamber centerline. The discrepancies 
observed in the results are attributed mainly to uncertainties in the initial spray 
size and velocity distributions and the droplet/wall impingement interaction, the 
single-step fast chemistry employed by the combustion model, and the deficiencies 
of the k - e turbulence model dealing with the strong streamline curvature. 

1.4 Dense Spray Models 

One of the important aspects in spray combustion modeling is the dense spray 
effects which include atomization process, drop breakup, droplet collision and co- 
alescence. Atomization process occurs on time and length scale too short to be 
resolved with practical computation grid sizes and time steps. Thus, atomization 
should be modeled as a sub-grid-scale process. To account for the dense spray ef- 
fects, the present study employs the drop collision & coalescence model[8] and the 
Taylor analogy breakup(TAB) model[9]. In the drop collision model, the probabil- 
ity distributions governing the number and outcomes of the collisions between two 
drops are sampled randomly in consistency with the stochastic particle tracking 
method. The TAB model utilizes an analogy between an oscillating and distort- 
ing droplet and a spring-mass system. The present breakup model is based on 
the reasonable assumption that atomization and drop breakup are indistinguish- 
able processes within a dense spray near the nozzle exit. Accordingly, atomiztion 
is prescribed by injecting drops which have a characteristic size equal to the nozzle 
exit diameter. Compared to Reitz’s model[34], the TAB model has several advan- 
tages in terms of no need to input the spray angle, an easy introduction of liquid 
viscosity effects, and explicit informations of distortion and oscillation effects on the 
interphase exchange rates of mass, momentum, and energy. For non-evaporating, 
evaporating, and burning dense spray cases, the predictions show a reasonably good 
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agreement with available experimental results in terms of spray penetration, drop 
size distributions, and overall characteristics of the evaporating and burning spray. 
Future studies include the detailed comparison with the local properties available 
in the experiment and the implementation of a volume-of fluid (VOF) formulation 
for resolving the liquid volume displacement effects. 


2. MATHEMATICAL MODELLING OF THE MULTI-PHASE FLOWS 


All the gas-phase and liquid-phase processes are modeled by a system of un- 
steady, two-dimensional(axi-symmetric) equations. The gas-phase equation is writ- 
ten in an Eulerian coordinate whereas the liquid-phase is presented in Lagrangian 
coordinates. The two-way coupling between the two phases is described by the 
interaction source terms which represent the rates of momentum, mass and heat 
transfer. These equations are given below. 

2.1 Basic Eulerian Equations 


2.1.1 Mean Flow Equations 

The density-weighted conservation equation of mass, momentum, and scalar 
variables in an Eulerian coordinate can be written as follows: 

+ d^ {pUi) = 5m,/ (21) 


dpui 

dt 


d 


dp 


+ ^-(pu.u,) = + 5 ., + s -,1 


dp<t> 

dt 




( 2 . 2 ) 

(2.3) 


where p is the time-mean density of the mixture, u, and u' x are the i component of 
the density- weighted mean and fluctuating part of the instantaneous velocity, <6 and 
<f>' are the density-weighted mean and fluctuating part of an instantaneous scalar 
quantities including the species concentrations and the total energy, p is the mean 
pressure, 5 $ and 5^,/ represents the gas-phase source terms and the interaction 
source terms due to the fuel spray, respectively. Detailed expressions for these 
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source terms can be found in Refs. [10, 11]. To close the system of equations, we 
need to model the unknown correlations, and u' t <p' . 


2.1.2 Turbulence Models 

The two-equation effective diffusivity model is used to represent the turbulent 
characteristics. In the eddy diffusivity models, the turbulent fluxes, u'u' and u '<£' , 
axe related to the mean flow gradients through the assumption of an isotropic eddy 
viscosity an d a constant turbulent Prandtl or Schmidt number. 


pu\u'j = 


. dui duj 2 . duk . 

Pt, d<f> 




(2.4) 

(2.5) 


The eddy viscosity (fit) appearing in (2.4) and (2.5) is defined in terms of a 
characteristic turbulence length scale(fc 3 ^/€) and a velocity scale (fc 1 ^ 2 ), so that fit 
is given by 

n, = C f p k j ( 2 . 6 ) 

The turbulent kinetic energy, k, and its dissipation rate, e, can be modeled from 
the turbulent transport equations: 


dpk d , ,, d / pt ^ dk — t j-dtti pt &P >7^ 

dt dxj 3 ’ dxj <Tk dxj 3 dxj p 2 dx } dxj 

dpe d , x d . . Pt^ de e . —f—fdui p t dp dp n e 2 . ^ 

-£-+-= — (pu,e) = -z—(p+ — )-z—+C(iT(P u i u ]sT~+-2A—A~)~ Cf2 P~u 
dt dxj 3 ' dxj a/dij fc 1 3 dx 3 p 2 dx } dx } k 

Here, terms involving jf- in (2.7) and (2.8) are inserted to account for variable- 
density effects[12]. These terms originally come from the pressure- velocity correla- 
tion in the Reynolds stress equation. For reacting flows, these terms should account 


partially for the expansion effect on the flow field due to heat release from combus- 
tion. 


2.1.3 Combustion Model 
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It is assumed that liquid fuel droplets act as distributed sources of fuel which 
evaporate to form a cloud of vapour. This implies that combustion process in spray 
flames can be treated as turbulent gaseous diffusion flames. Experimental evidence 
for this assumption can be found in Refs. 13. An idealized approach for physically- 
controlled diffusion flames is to invoke a fast-chemistry assumption which the chem- 
istry is sufficiently fast and intermediate species do not play a significant role. In 
the turbulent diffusion flame model, the influence of turbulence on combustion is 
taken into account by relating the fluctuations of mass fractions. This implies that 
fuel and oxidizer can coexist in the same place but at a different time. The most 
convenient way to include the effect of turbulent eddies on thermochemical proper- 
ties is via the introduction of the probability density function(pdf), P(£, x,). This 
function contains information of both mean(f) and variance of ( g = (/ — /) 2 ) of the 
mixture fraction. These variables f and g can be obtained by solving the transport 
equations. The density-weighted mean values(<£) of any property are evaluated by 
convoluting the property functions with a probability density fuction, P(£, £,): 

3=fn()P&*iW ( 2 . 9 ) 

Numerous probability density functions are available in the literatures. The present 
study adopts the /? pdf which is known as the widely applicable one [11,12]. The 
detailed numerics of the 0 pdf are well described in Appendix B. The double - 
delta pdf procedures are also implemented in the program. 

A modified eddy breakup model[14] is optionally incorporated in the present 
computer code. Using this model, the reaction rate is determined as follows: In an 
irreversible sigle-step chemical reaction, the mixing-controlled reaction rate[14] is 
given by 


Rmix = Ap^min(Y f , y) (2.10) 

where A is a model constant; s is the stoichiometric oxidant/fuel ratio; Yf and 
Y 0 are the mass fraction of the fuel and the oxidizer. To account for the ignition 
delay time, the chemical kinetics need to be considered. The chemically controlled 
reaction rate, P c /, e , is given by the usual Arrhenius formula[15]. 
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( 2 . 11 ) 


Rcht 



w 0 } 


The reaction rate, R/ u is determined from either of the mixing rates of the reactants 
or the chemical reaction rate, whichever slower. 


Rf u = min(R mix > R ch <) (2-12) 

For simple one-step reaction of the hydrocarbon-air mixtures, there are five species 
participating the mixture composition. Once the mass fractions of fuel and oxdizer 
have been determined, the remaining species can be easily determined from the 
stoichiometric relations described in Appendix C. 


2.2 Basic Lagrangian Equations 


'W 


2.2.1 Droplet-Phase Equations 


In this study, the spray is described by a discrete particle method formulated 
on a Lagrangian frame. This is essentially a statistical approach and requires track- 
ing a sufficiently large number of computational particles. Each computational 
particle represents a number of droplets having equal location, velocity, size, and 
temperature. The governing equations for these are . 


doc | 



(2.13) 


dpdVi _ Uj + - V i p J 

dt Ti 

dr p TT 7. ev 

dt 4nr p 2 pd 


(2.14) 

(2.15) 


and 


dTd 

dt 


Ql 

TTljiCp^d 


(2.16) 


[n equation (2.14), Fb t represents the body force terms such as the gravity force 
and the centrifugal force. The particle relaxation time r, can expresses as: 

r" 1 = | — Co\u t + u t ' - u;| (2.17) 

o r t> 
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Cd is the drag coefficient given by 


C D = P-( l + s«e, } ); ^ Re, <1000 

Afijj y 

and 0.424; for Re p > 1000 


In which 


Re p = 


I Ui +u,' - v t \pd t 


(2.18) 


(2.19) 


In equation (2.15), the droplet evaporation rate is given by the Frossling correlation 
[16,17] 

m ev = 2xd p (pD)(l + O.ZRe p ^ScJ)ln(l + B m ) (2.20) 

In equation (2.16), the droplet temperature, which is assumed to be constant within 
the droplet, is found by using the heat energy Ql, • 


Ql — 47rr J> 2 <3c - rn ev L 


( 2 . 21 ) 


where L is the latent heat of vaporization, and Q c is the heat conduction rate to 
the droplet surface per unit area. Q c is given by the Ranz-Marshall correlation 

= 2R(T z 7 k ) (1 + 03flef W) foO^) 
d p 

The Schmidt number, Prandtl number and mass transfer number are defined re- 
spectively as _ 

c _ JL- Pr -0 £l 

SCj — _ i rr 


Y„ —YqQ y _ Pj_ 

Bm ~ 1 -Y, ’ P 


(2.23) 


The values of thermodynamical properties of gas such as K , C p , D etc. are 
highly dependent on the temperature and fuel vapor mass fraction at which they are 
evaluated. A "one-third rule” [18] that utilizes a reference temperature equal to the 
droplet surface temperature plus one-third of the difference between the surrounding 
gas and droplet surface temperature is used. The same procedure is applied to the 
reference value for the fuel vapor mass fraction, in which Y, is obtained from 


p w'- 1 


n = [l + (^--l)irrl 


(2.24) 
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Here Y, and P v axe the mass fraction and the fuel vapor pressure at the droplet 
surface, and W f and W a are the molecular weights of fuel and mixture, respectively. 
For a given Td , Pv is estimated from the JANAF data bank [19]. The two-phase 
interaction source terms in the gas-phase governing equations are described in Ap- 
pendix A. The droplet distribution models for the dilute sprays are also described 

in Appendix C. 

In case of the droplet passage through the plane of symmetry, the droplet 
with same instantaneous properties and physical dimensions, but the mirror image 
velocity vector, is injected into the flowfield. On impingement on a wall, the droplets 
are assumed to bounce back with the reduced momentum(lO). 

2.2.2 Turbulence/Droplet Interactions 

In this study, the spray is described by a discrete particle method formulated 
on a Lagrangian frame. To account for turbulence dispersion, we follow the concept 
of [7] of combining a normal (Gaussian) probability distribution for each computa- 
tional particle. The instantaneous location of each computational paiticle is calcu- 
lated by a stochastic Lagrangian tracking scheme. The governing equation for each 
computational particle is 


dv k u k -Vk , r 

+ Fbk 

dt r k 

(2.25) 

*11" 

II 

<2 

(2.26) 

-> = l£-Cp.\ Uk - Vk \ 

* 8 p f d„ 

(2.27) 


The location calculated by the above equations only represents the mean of each 
particle’s corresponding probability function. The variance of each parcel pdf has 
to be calculated and the combined pdfs then represent the statistical distribution 
of particles with turbulent dispersion effects. To estimate the variance of the parcel 
pdf due to the turbulent particle dispersion, the turbulence-induced displacement 
and velocity can be splitted from equations (2.25) and (2.26): 

dv\ u k ~ k (2 OS) 

dt ~ r k 
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(2.29) 


dx' 

dt 


k t 

= V k 


With the isotropic turbulence assumption, each component of u' k is randomly 
chosen from a Gaussian distribution with standard deviation u'* rmj = \J\k- We 
first choose A t kt as the time step of the i th interaction within the k th eddy, which is 
smaller than the eddy lifetime, and integrate equations (2.28) and (2.29) to update 
particle fluctuating locations and velocities. 


■At| 


x' kt = u' krma A t ki + (*'*«_!) - )(1 - e ) (2.30) 

v'ki = u'krms + i v> k(i- 1 ) - “VwjK* 1 " 1 ’ (2.31) 

We then sum up the m steps for which the particle fully interact with the k th eddy, 


53 i = A t k 


(2.32) 


i=i 


The change of variance of a computational particle pdf within the k th eddy is 
represented by a characteristic mean squared dispersion in the form: 

m 2 

V = **-■’ + (£*'*) (2-33) 


1=1 


In equation (2.33),<7fc_i is the existing variance of the particle pdf at the begining of 
the interaction within the k th eddy. Since the time step within each turbluent eddy 
is fixed, the number of interaction within the eddy, m, varies across the calculation 
domain, the choice of time step A t kl and the related issues are discussed in detail 
in [24]. Figure 2.1 well describes this eddy interaction with the particles. The 
present procedure is easy to program and requires less computer memory. For 
each computational particle, we just need to store x' k i, u' krms ,v' kt , and cr k 2 . This 
procedure when implemented in the current time-marching numerical method is 
somewhat different from the method of [7] in which the calculation of the current 
variance of each particle pdf is summed over the entire history of the effective 
time constants. In their recent study, truncation of unnecessary time history terms 
and the associated errors was discussed and additional computational efficency was 
obtained [35]. 
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When convoluting pdf for a group of computational particles, the variances 
of eq.(2.33) must be normalized according to the total number of paxticles. The 
normalized particle variance can be written as 

», k = Kjji- (2.34) 

Here, represents the statistical uncertainty in the mean particle position, K is 
the correction factor to account for undersampling, and Nt is the total number of 
computational particles. When symmetry and reflective boundary condition exist 
in the calculation domain, a cumulative pdf distribution at any point in coordinate 
y, y is the distance from the particle to the axis or the reflective boundary, may be 
defined as: 3 

P(y)= f -=U-e ""U dy (2.35) 

J-y V2tT<7 y k 

Here, y p is the instantaneous location of computational particles. After integration, 
the symmetric cumulative distribution function takes the form, 

P(») = 0.5(er/(^^)+er/(^±^)] (2.36) 

VtVyk VtOyk 

where 



In accordance with the approach of Litchford and Jeng[7], when the mean po- 
sitions of computational particles is calculated by the deterministic tracking(u* in 
Eqs.(2.25)-(2.27) is the mean gas velocity), this approach is described as the deter- 
ministic dispersion width transport(DDWT) model. For tracking using stochastic 
sampling of gas-phase turbulent velocity fluctuations^* in Eqs.(2.25)-(2.27) is the 
instantaneous gas velocity), the approach is described as a stochastic dispersion 
width transport(SDWT) model. 

In the point delta function SSF model, the turbulence effects on droplet dis- 
persion are simulated by a Monte Carlo method in the sense that a fluctuating 
velocity u'*, where each component of u ' * is randomly chosen from a Gaussian dis- 
tribution with standard deviation J~\k, is added to the mean gas velocity. Thus 
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the turbulence is assumed to be isotropic. This type of simulation for the turbulent 
dispersion of droplets has been extensively used previously [20,21,22] for statisti- 
cally stationary turbulent dispersed flows. Main differences in the implementations 
are the methods used to specify turbulence eddy properties and the methods for 
choosing the time of interaction of a particle with a particular eddy. The details of 
simulation procedures and also of various aspects associated within the interaction 
times can be found in Ref. [11,38]. 

2.2,3 Drop Breakup and Collision 

The present study employs the TAB (Taylor Analog}' Breakup) model pro- 
posed by O’Rourke and Amsden[9]. This model is based on an analogy between an 
oscillating and distorting droplet and a spring-mass system. The restoring force of 
the spring is analogous to the suface tension forces. The external force on the mass 
is analogous to the gas aerodynamic force. The damping forces due to liquid viscos- 
ity are introduced to this analogy. Compared to Reitz’s model[34], the TAB model 
has several advantages in terms of no need to input the spray angle, an easy intro- 
duction of liquid viscosity effects, and the explicit informations of distortion and 
oscillation effects on the interphase exchange rates of mass, momentum, and energy. 
The major limitation of the TAB method is that only one oscillation mode can be 
tracked. However, in reality there exist many such modes in the Taylor analogy. De- 
spite this limitation, good agreement between numerical results and experimentally 
observed bag/stripping breakup times has been reported. The droplet oscillation h 
breakup calculations require two normalized particle array s( deformat ion and oscil- 
lation) which can be determined by the equation for the acceleration of the droplet 
distortion parameter. Occurance of droplet breakup, the Sauter mean radius(SMR), 
and oscillation velocity for the product drop depend on these two parameters and 
Weber number. The radius of the product drops is then chosen randomly from a 
chi-squared distribution with calculated SMR. Following breakup, the product drop 
has the same temperature with the parent drop, and its deformation and oscillating 
parameters are set to zero. 

The drop collision model suggested by 0'Rourke[8] is employed to calculate 
collision and coalescence among the dispersed liquid phase. The collision routine is 


operatinging for the pair of particles if, and only if, they are in the same computa- 
tional cell. For the collision calculation, the drops associated with each computation 
parcel axe considered to be uniformly distributed throughout the computational cell 
where they are located. For all parcels in each computational cell, a collision fre- 
quency between drops between the paice\{parceli) of larger drop radius^) and the 
parcel(parce/ 2 ) of smaller drop radius(r 2 ) is obtained from the relationship in terms 
of the number of drops in parcel 2 , the relative velocity between parcel i and parcel 2 , 
the area based on n + r 2 , and the volume of computational cell. The probability 
with n collisions is assumed to follow a Poisson distribution based on a collision 
frequency and the computational time step. Using the probability informations, 
the collision impact parameters are stochastically calculated. If the collision impact 
parameter is less than a critical impact parameter, the outcome of every collision 
is coalescence. In opposite case, each collision is a grazing collision. The critical 
impact parameter depends on the drop radii, the relative velocity between drops, 
and the liquid surface tension coefficient. 

3. NUMERICAL MODEL AND BOUNDARY CONDITIONS 

The present method is based on the operator splitting technique[23] attempting 
to reach accurate transient solution after prescribed predictor- corrector steps for 
each time-marching step. The previous multiple pressure-correction method[11.27) 
is extended to to handle the strong nonlinear couplings arising in the multi-phase, 
fast-transient, and reacting flows. This method is non-iterative and applicable to 
all-speed flows. The additional scalar conservation equations such as species, and 
enery are incorporated into the same predictor-corrector sequence. Discretization of 
the gas phase governing equation uses the finite volume approach. To enhance the 
numerical stability, the implicit Euler scheme is employed in differencing the tempo- 
ral domain. All the dependent and independent variables are stored at the same grid 
location and the variables at the finite control volume boundaries are interpolated 
between adjacent grid points. The discretizations have been performed on a general 
non-orthogonal curvilinear coordinate system with a second order upwind scheme 
for convection terms and the central differencing scheme for diffusion terms. The 
resulting discretized equations were solved by a conjugate gradient (CGS) solver. In 
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the present algorithm, each time step is divided into a one-predictor/two-corrector 
sequence. The strong coupling terms between particle and gas are evaluated by the 
same time splitting technique. Implicit coupling procedures are used to treat mo- 
mentum exchanges to avoid the small timesteps. The unsteady solution procedure 
described above is somewhat different from the conventional PSIC(particle source 
in cell) procedure[25] in which global iterations are required. The method used here 
is non-iterative and time-accurate. 

To improve the convergence and the numerical stability for the fast transient 
spray-combusting flows, all transport equations except the continuity equation axe 
expressed in the advective form. 


pn+l -p n d , ,n+l _ <j 


(3.1) 


l A t* +V,) dxi 




dxi 9 dx{ 




(3.2) 


By using operator-splitting method, the transport equations with the predictor- 
corrector procedure can be discretized as follows: 


(af Predictor step 
Momentum(u*): 


p n u? 


{ At _ Ao)< = Hn{U ' ] _ A ' pn + + + a/ 


(3.3) 


Scalar(^*): 



-BoW = J n (4>') + s; + s; tt + 


p n <? n 

At 


(3-4) 


Here, the operators A 0 , B 0 , H, and J are constructed from the third-order upwind 
scheme for the convection terms and the central differencing scheme for the diffusion 
terms. To improve the numerical stability in multi-phase reacting flows, A 0 and B 0 
may include the coefficient matrix resulting from the implicit treatment of the strong 
nonlinear source terms such as chemical reaction rates, turbulence source terms and 
multi-phase interaction source terms. The quantities S" ( , S u r |, 5^ , and t are 
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determined from the existing flow fields. The general scalar dependent variables, (f> 
may represent the enery, the mass fraction, and the turbulent transport quantities. 
In this stage, the velocity field(u*) does not satisfy the continuity equation. The 
temperature T * is calculated from the flowfield(energy, species, momentum) at the 
predictor step. 

(b) First Corrector step 


Moment um(u**): 

A new flowfield(u**, p*, p *) are sought to satisfy the continuity equation: 

- p") + A^-u") = SZ., 


(3.5) 


and the discretized momentum equations axe: 

( ^ = **(«*) - + s l + s l ,, + 


(3.6) 


Continuity equation (3.5) can be rewritten as: 


- p") + A.(^(u” - uj)l + A,[(p- - ,>>"] = -A.(/u-) + 5" ,, (3.7) 


Subtracting Eq.(3.3) from Eq.(3.6) gives the velocity correction equation. 


p’W - «•) = -p"Du"\A,(p' -P")l 


(3.8) 


Here, 

Du' = (£ t -Al)" 

Equations (3.7) and (3.8) axe now used to derive the pressure correction equation. 
Thus, talcing the divergence of Eq.(3.8) and substituting into Eq.(3.7) yields 




1 ^ + A i (^)-A i (/>"r»u"A i )](p--p") = 


AtRT* ' RT' 

-A ,(A*) + 5- ,, + (^ - ~)^- t + A.[( ^ - ^)P n «\ (3.9) 


1 1 ,p" . 1 1 


Equation (3.9) can be solved for the corrected pressure, p*. The density(p’) is cal- 
culated from the equation of state. The velocities(u**) axe computed from Eq.(3.8). 
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Scalar( <£**): 

These new flowfield(u’V) satisfying the continuity equation (3.5) are used to 
update the B coefficient. 


it 

P <P 


( - BDir = •/'«•) + s; + s;, + 


(3.10) 


The temperature T" is calculated from the corrected flowfieldfenegy, 

species, momentum). 

fc) Second Corrector step 

Momentum(u*** ): 

A new flowfield(tt***, p**, p**) are sought to satisfy the continuity equation: 


1.( P " - p -) + A,(p"u‘") = SI,, 
Subtracting Eq.(3.5) from Eq.(3.11) yields 

£(/>" -/.*) + A, [/(a*" - «**)] + A,[(p" - P>*“] = 0 

* * V 


(3.11) 


(3.12) 


The discretized momentum equations are: 

( £ - K)<" = **(«” ) - A.p" + S" + 5" ,, + 
Subtracting Eq.(3.6) from Eq.(3.13) yields 


At 


(3.13) 


*** . +* 

Ui - u, 


= Du-[(A : - A n 0 )u •• + JT(u " ) - H n (u*)] - Du'[A t (p” - p)\ (3.14) 


Here, 


-1 




By substituting Eq.(3.14) into Eq.(3.12), the corrected pressure equation is 
obtained: 


[_J_ + 4( |L ) . a ,( p - D „-A,, r - / ) = 

&,{/)• Du‘[( Al - AJJa" + H'( u”) - ff"«)]] 
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(3.15) 
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RT* 


1 

RT 


^s +A - [( 
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RT* 


1 

RT 


)p*u*] 


Using the corrected pressure(p**), the density(p**) is obtained from the equation 
of state. The velocities(u***) are computed from Eq.(3.15). 


Scalar^***): 

These new flowfield(u" V) satisfying the continuity equation (3.11) are used to 
update the B coefficient. 


_ 71 

(~r B 

v At 


W = r m (<r) + s;* + s; tl + 


in 

P <P 
At 


(3.16) 


The temperature T*** is calculated from the updated flowfield(energy, 
species, momentum). The updated fiowfield (/>**, p**,u***, T***, and <j>***) are taken 
to represent the field values at the next time step(n+l). This completes the se- 
quence in the solution of the equation over the time-step. 

In the present pressure-based method, the left hand side of the corrected pres- 
sure equations written in Eq.(3.9) and Eq.(3.15) involve a convection term which 
can properly takes into account the hyperbolic nature of supersonic flows, and en- 
ables capturing shock waves. A recently developed non-iterative PISO method[39] 
does not include this convective term due to the inconsistent treatment of density 
in the momentum equations. Compare to our previous pressure-based method[27], 
this new method allows the consistent discretization of continuity equation and 
is more suitable for the fast-transient reacting flow at all speeds. For the steady 
state calculations, the present procedure can be simplified by freezing the coefficient 
matrix(A 0 , B 0 ). 

For the subsonic inlet boundary, the entropy and the total pressure are spec- 
ified. The axial velocity components are obtained by the extrapolation and the 
vertical velocity components are determined by enforcing the vorticity to vanish at 
the upstream boundary. At symmetry, the normal grdients of all scalar variables 
and the radial velocity component are zero. At the supersonic outlet, all dependent 
variables are extrapolated from the interior. The wall was assumed to be adiabatic. 


4. RESULTS AND DISCUSSIONS 
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To evaluate the present dispersion width transport model and to calibrate 
the stochastic simulation of particle-turbulence interactions, the computations were 
performed for the solid particle dispersion in nearly- homogeneous turbulence and a 
particle laden round jet in inhomogeneous turbulence. The validation cases for the 
dense spray models include non-evaporating, evaporating, and burning sprays. 

4.1 Turbulent Particle Dispersion 

4.1.1 Nearly- Homogeneous Turbulent Dispersion 

The particle dispersion experimental setup of Snyder and Lumley [28] in a 
grid-generated turbulent flow was used for evaluating the present dispersion PDF 
transport model. Particle densities and sizes are chosen to examine the phenomena 
in which the eddy lifetime controls interaction times (46.5 /im diameter hollow- 
grass) , the transit time controls interaction times (87.0 //m corn pollen), or the 
controlling-interaction times undergo transition from transit time to eddy life time 
(87.0 /jm solid glass). In this experiment, fluid turbulence intensities and length 
scale information were measured. The particle calculations were started at the ex- 
perimental particle injection point of x/m = 20 (m is a 2.54-cm-square mesh). The 
particle velocity was assumed equal to the mean fluid velocity of 6.55 m/sec. For 
the delta function SSF computations, 5,000 computational particles were sampled 
to calculate the resulting mean squared dispersion with respect to time. For the 
DDWT computations, a single parcel in a deterministic trajectory along the center- 
line was sampled to evaluate the mean squared dispersion representing the variance 
of the parcel PDF by using the related parameters for each eddy interaction. 

Figure 4.1 shows comparison of the predicted and measured particle dispersion 
with respect to time. The DDWT results show good agreement with the SSF results 
for light, medium, and heavy particles. Both models also show favourable agreement 
with the experimental data. These numerical results indicate that the DDWT 
model with a single computational parcel following the deterministic trajectory 
demonstrates the efficiency, the accuracy, and the overall prediction capability for 
this nearly-homogeneous turbulent flow. 

4.1.2 Inhomogeneous Turbulent Dispersion 
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The next example problem is a paxticle laden round jet [29] in which the tur- 
bulence is inherently inhomogeneous. The turbulent gas-phase transport properties 
are provided by using the k-e model. Figures 4.2 and 4.3 show the particle concen- 
tration profiles of the delta function SSF model and the SDWT model for 50 and 
200 computational parcels, at various levels of the correction factor, and at several 
axial locations. 10,000 particles are sampled for the delta function SSF computa- 
tions. Using the 10,000 particles in the SSF model, there is still evidence of slight 
undersampling. However the distribution is relatively smooth and is taken here as a 
good approximation to the theoretical profile. The 50 parcel case of SDWT model 
shown in Figure 4.2 is very sensitive to the level of the correction factor especially 
for upstream regions due to undersampling. By increasing the correction factor, K 
in eq.(13), the uncertainty level in the mean increase the dispersion and smooth the 
profile considerably. In Figure 4.3, the zero correction factor case(K=0) corresponds 
to the delta function SSF case using 200 computational particles. The computed 
profile of the SSF case for 200 particle samples is very irregular and shows oscilla- 
tory distribution. The 200 parcel case of SDWT model shown in Figure 4.3 is less 
sensitive to the correction factor since there is less uncertainty in the mean because 
of increased s am pling. In Figure 4.4, the SDWT results with 200 parcels and K=4 
shows favourable agreement with the delta function SSF with 10,000 computational 
particles. These numerical results clearly indicated that the SDWT model has the 
capability of accuratly representing dispersion in inhomogeneous turbulent flows 
with improved efficiency over the delta function SSF model. 

4.2 Dilute Spray Combusting Flows 

The present numerical model for the multi-phase turbulent reacting flows has 
been tested by applying it to predict the local flow properties in two axisymmetric, 
confined, swirling spray-combusting flows[36,37]. 

4.2.1 Hollow-Cone Kerosene Spray Flames 

The combustor geometry of the first test case is shown in Fig. 4.5. Experimen- 
tal data for temperature, axial and tangential velocity components were obtained 
from measurement of Khalil et. al. [36] . The inlet conditions and the initial droplet 
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size distribution axe given in Table 1. Liquid kerosene was used as fuel and the 
air/fuel mass ratio was fixed at 20.17. 

In the present study, two swirling numbers(S=0.72 and 1.98) were considered 
to investigate the influence of swirl on the droplet evaporation & burning charac- 
teristics. Fig. 4.6-4. 8 show the general flow pattern such as the predicted droplet 
trajectories, velocity vectors, and temperature contours of two swirl cases. In the 
lower swirl case(S=0.72), large portion of droplets survive in the central recircu- 
lation zone and continue to evaporate in the fax downstream region. In the high 
swirl case(S=1.98), most of small droplets axe trapped in the recirculation zone and 
evaporate there, producing intensive burning and high temperature in this region. 
With increasing swirl, the droplet spreading increases due to the droplet dispersion 
and the increased particle centrifugal force term. In addition, the larger central re- 
circulation zone corresponding to the higher inlet swirl is contributed to recirculate 
more hot combustion gas from downstream and to increase the temperature at near 

inlet regions. 

The predicted and measured temperature profiles for two swirl cases axe shown 
in Fig. 4 9 and 4.10. The siginificant discrepancies in near- wall regions mainly 
result from the uncertainties of droplet/wall impingement process. However, the 
deviations in other locations are associated with the deficiencies of turbulence and 
combustion models, the unreliable informations of the inlet droplet size & veloc- 
ity distribution, and the potential errors in inlet swirl profiles and inlet turbulence 
length scale. It is observed that the temperature profiles of the high swirl case 
axe more uniform than those of the low swirl case. Radial profiles of axial velocity 
velocity for S=0.72 and 1.98 axe shown in Fig. 4.11 and 4.12. The present numer- 
ical model underpredicts the magnitude of the reverse flow velocities. The poor 
performance of the present numerical model in predicting the size of central recir- 
culation zone and the reverse velocity is partly attributed to the deficiency of k - e 
model based on the isotropic assumption. The predicted and measured tangential 
velocities for two swirl cases are presented in Fig. 4.13 and 4.14. The significant 
deviation close to the inlet is likely caused by the incorrect distribution of inlet 
swirl velocities. In the present study, the inlet swirl velocities are obtained from the 


estimated axial angular momentum flux. The rapid decay of the tangential velocity 
to the solid body rotation close to the centerline could be tied with the errors in 
the prediction of reverse velocities. 

4.2.2 Spray Flame with a Rotating-Cup Atomizer 

Fig. 4.15 shows the liquid-fueled combustor with a rotating cup atom- 
izer[37] which is capable of producing a near-monodisperse spray. This small near- 
monosized spray allows the relatively accurate representation of the droplet initial 
conditions and eliminates the uncertainties in the droplet/ wall impingement process. 
The gas-phase inlet boundary conditions and the droplet intial con ditions( droplet 
size k velocity distributions) were estimated from measurement of El-Banhawy and 
Whitelaw[37]. Liquid kerosene was used as fuel and the fuel/air mass ratio was 
fixed at 0.0228. In present study, computaions were carried out for a test condi- 
tion with the swirl number(S=1.2) and the droplet diameter(r* = 47/jm). Since 
the experimental data provided the limited informations for gas-phase inlet and 
droplet injector conditions, there are some potential errors associated with the in- 
let swirl profiles, and the initial distributions for droplet size and velocity. In this 
test case with small droplet size, considerable change in droplet velocity can occur 
in the short distance between the injection location and the measurement station; 
therefore, some error might be introduced in the droplet initial velocity specifica- 
tions. Based on the experimental data, the droplet size distributions were specified 
with the near- monosized droplets( 47/im, 70 % fuel mass flow rate) and the smaller 
satellite droplets (24 /im, 30 % fuel mass flow rate). 

The predicted droplet trajectories, velocity vectors, and temperature conturs are 
shown in Figure 4.16-4.18. It is observed that two high-temperature regions exist 
in the shear layer around the recirculation zone, and the main flame region around 
the central recirculation zone and downstream of the fuel spray. Numerical results 
indicate that these high-temperature regions are characterized by the trapping of 
smaller droplets, high evaporation rate, and intensive turbulent mixing and chem- 
ical reaction. Figure 4.19-4.21 show the predicted and measured radial profiles 
of temperature, mass fractions of carbon dioxide ( CO 2 ) and oxygen(C> 2 ) at four 
axial locations. The numerical results shows the qualitative agreement with the 


experimental data. However, quantitative differences exist in the profiles of tem- 
perature and corresponding chemical concentrations. The sigificant discrepancies 
close to the inlet(X/D=0.254, 0.510) are mainly attributed to the incorrect speci- 
fication of swirl velocity profile, the insufficient turbulent mixing predicted by the 
k - e turbulence model, and the neglect of intermediate species such as soot and 
carbon monoxide(CO) associated with the single-step fast chemistry model. In 
these regions, the overpredicted CO 2 mass fractions can be partly explained by the 
existence of the high CO mass fractions^ 10%) observed in measurement. Due 
to the relatively high heat release rate of C0 2 , the calculated temperatures are 
much higher than the measured values. At the far downstream regions(X/D=2.9) 
of the spray flame, the differences between the calculated and measured temperature 

values decreases to around 100 K. 

4.3 Dense Spray 

4.3.1 Non-Evaporating Solid-Cone Spray 

The solid-cone spray measurements of Hiroyasu and Kadato[30] were used to 
validate the present numerical dense spray model which includes collison, coales- 
cence, and breakup models described above. Liquid fuel is injected through a single 
hole nozzle into constant pressure, room-temperature nitrogen. Spray tip penetra- 
tion and drop sizes were measured from photographs of the backlighted spray. The 
test conditions are given in Table 2 (SMD is the average over the spray cross-section 
65 mm downstream of the nozzle). The nozzle diameter was 0.3 mm and the present 
computations used tetradecane for the liquid fuel(the experiments used a diesel fuel 
oil with physical properties close to tetradecane). 

A computational domain of 20 mm in radius and 120 mm in length was dis- 
cretized by a 25 radial and 45 axial grid. The mesh spacing was nonuniform with 
refinement on the centerline and close to the injector. The smallest cell is 0.5 mm 
radially and 1.5 mm axially. Since this dense spray calculation is sensitive to the 
grid resolution, the fine grid was used to obtain a grid-independent solution. The 
number of computational parcels at steady-state conditions was between 1000 and 
1500, and the number was varied with the back pressure. The present numerical 
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results did not change appreciably when this parcel number was varied. The initial 
turbulent quantities were assumed as the small values(fc = 1 x 10 3 m 2 /s 2 ,e = 4 x 
10 — *m 2 /s 3 ). The numerical results were insensitive to these initial values. 

The spray parcel distribution for three sprays is shown in Figure 4.22. This 
plot indicates that the spray tip penetration and the core length decrease with the 
increase of the gas density. Figure 4.23 shows the predicted and measured spray tip 
penetration versus time. It can be seen that there is reasonably good agreement 
between the prediction and the measurement . In the present computations, the 
spray tip was defined to be the location of the leading spray drop parcel. It is 
necessary to note that a far-field spray penetration is not a sensitive indicator of 
model performance. Previous studies[26,31] indicated that a far-field spray penetra- 
tion is mostly influenced by the turbulence diffusivity. However, a near-field spray 
penetration could be more sensitive to the physical submodels such as breakup and 
collision. Figure 4.24 shows the variation of SMD with axial distance from the 
injector. The three solid data at 65 mm correspond to the measurements. The 
computed drop size is a time average over the spray cross-section at each axial lo- 
cation. At the nozzle exit, the drop diameter is equal to the nozzle diameter, 0.3 
mm Generally these curves can be broken into two sections. Close to the injector, 
the drop size decreases rapidly due to drop breakup. Further downstream, the drop 
size increases gradually due to drop coalescence. In the low gas pressure case(l.l 
MPa), the drop size remains relatively uniform after initial breakup region and then 
increases slightly in the far-downstream region. For the high pressure cases(3.0 and 
5.0 MPa), the drop size increases largely in far-downstream region, because higher 
gas densities promote collisions and coalescence. This trend is also observed in the 
measuments. The predicted drop sizes at 65 mm are qualitatively agreed with the 
experimental data for all three cases. The discrepancy could be associated to the 
fact that the experimental sprays were pulsed while the computations assumed a 
constant pressure injection for the entire computational time period. 

4.3.2 Non-Evaporating Hollow-Cone Spray 

The hollow-cone spray tip penetration data of Shearer and Groff|32] have been 
used for the model validation. In the experiment, the liquid is injected into quiescent 
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room-temperature nitrogen at P = 550 kPa. The numerical timestep is used as 
2.5^3 and about 2000 spray parcels are used in the computation. The experimental 
spray cone angle is 60 degrees, and the flow rate 0.0165 mL/injection with four 
pulses, each of duration about 0.58 ms. The computational injection velocity are 
set to the experimental spray tip velocity(60m/s) measured from the movie pictures 
in the early stage of the injection. The test condition is listed in Table 3. 

Figure 2.25 shows the spray parcel distribution and the velcocity vectors, and 
the predicted and measured spray tip penetration versus time. The numerical results 
indicate that turbulence has a relatively small effect on penetration in a hollow-cone 
spray because radial spreading due to inertia is dominant. The gas velocity vectors 
indicate the presence of a vortex near the head of the spray, which curls the spray 
tip toward the outside of spray. A substantial region of strong inward flow in the 
center of the cone near the injector was also observed. These flow patterns and spray 
shapes compared quite favorably with the experimental observations[32]. In Figure 
2.26, the predictions reasonably agree with the experimental spray tip penetration. 

4.3.3 Evaporating and Burning Solid-Cone Spray 

The evaporating and buring solid-cone spray measurement of Yokoda et. al. [33] 
have been used to validate the present numerical dense spray model. Liquid 
fuel( tridecane) is injected through a sigle hole nozzle into high-pressure, high- 
temperature nitrogen or air. The test conditions for evaporating and burning sprays 
are given in Table 4. The nozzle diameter was 0.16 mm. A computational domain 
of 20 mm i n radius and 100 mm in length was discretized by a 21 radial and 44 
axial grid. The mesh spacing was nonuniform with refinement on the centerline and 
close to the injector. The number of computational parcels at steady-state condi- 
tions was between 500 and 700. Due to the numerical reasons, the initial turbulent 
quantities were ass um ed as the small values. The upstream boundary is treated as 
a solid wall, and other boundary are treated as open boundaries. 

Figure 2.27 shows the spray parcel distribution and the contours of the fuel 
mass fraction for evaporating sprays. These results show that the spray penetra- 
tion increase with respect to time at early period of injection, however the pene- 
tration become nearly constant after t — 0.2 ms due to evaporation. Even though 


the liquid drop does not penetrate more, the evaporated fuel vapor continuously 
penetrate with respect to time. Comparisons of the computed and experimental 
spray penetration versus time are shown in Figure 2.28. The present spray penetra- 
tion distance agrees well with the measured results[33]. Figure 2.29 and 2.30 shows 
the spray parcel distribution, the contours of the fuel mass fraction, temperature, 
and oxygen mass fraction at different times of injection for burning sprays. The 
computed configuration of a burning spray flame has the overall agreement with 
the measure ones. In the experimental study, a considerable level of soot was ob- 
served near the spray tip where the equivalence ratio is low and the temperature 
is high due to the progressed turbulent mixing. Therefore, the soot model should 
be incorporated to improve the prediction capability of the present burning dense 
spray model. Future studies may include the detailed comparison with the local 
properties available in the experiment. 

5. CONCLUSIONS 

The numerical models have been developed for the analysis of dilute and dense 
spray-combusting flows. From the present numerical studies, the following conclu- 
sions are drawn in general: 

1. Present implementation of the dispersion width transport model has success- 
fully demonstrated the capability of accuratly representing dispersion in nearly- 
homogeneous and inhomogeneous turbulent flows with improved efficiency over 
the delta function SSF model. 

2. A numerical model for the prediction of the statistically stationary spray- 
combusting flows is evaluated by comparison with the available experimental 
data. The present numerical procedure correctly predicts the general features 
of spray-combustion flows and yields the qualitative agreement with experi- 
mental data. However, quantitative differences exist especially at near-burner 
locations, at near- wall regions, and along the combustion chamber centerline. 
The discrepancies observed in the results are attributed mainly to uncertainties 
in the initial spray size and velocity distributions and the droplet/ wall impinge- 
ment interaction, the single-step fast chemistry employed by the combustion 
model, and the deficiencies of the k — 6 turbulence model dealing with the strong 
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streamline curvature. To improve the prediction capabilities of the present 
numerical procedure, the future works must include the consistent studies of 
non-evaporating, evaporating, and burning sprays by utilizing the non-isotropic 
turbulence model such as the algebraic stress model and the second-moment 
closures, and the multi-step finite chemistry model. 

3. For non-evaporating, evaporating, and burning dense spray cases, the predic- 
tions show a reasonably good agreement with available experimental results in 
terms of spray penetration, drop sizes, and overall configuration of a burning- 
spray flame. To improve the prediction capabilities and efficiencies of the nu 
merical and physical models, future works must include the extensions of the 
dispersion width transport model to non-evaporating, evaporating, and burn- 
ing dense sprays, the incorporation of supercritical vaporization model, the in- 
corporation of soot model and further refinement of atomization and breakup 
models. 


6. RECOMMENDATIONS 

The following recommendations are intended as suggestions for improvements 
and extensions of the present spray combustion modelling. The numerical and 
physical modelling studies are need for the following tasks. 

• Implementation of computationally efficient parcel PDF approach for multi- 
phase, turbulent, evaporating, and combusting flows. 

• Development of strong interphase coupling procedure by combining multiple 
pressure correction procedure and Volume of Fluid(VOF) method. 

• Im plementation of equilibrium and non-equilibrium chemitry packages for effi- 
cient transient reacting flow calculations. 

• Optimization and adaptation of breakup and coalescence procedure. 

• Atomization modeling in conjunction with multi- step pressure correction 
methodology. 

• Incorporation of turbulence modulation effects by droplets. 

• Incorporation of wall/droplet impingement process. 
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• Incorporation of supercritical vaporization model. 

The validations and the applications of the proposed spray combustion models will 
be consistently studied for the following cases: 

• Benchmark solution for non-evaporating, evaporating, and burning dense 
sprays. 

• Unsteady flame propagation in a two-dimensional spray with transient droplet 
vaporization. 

• Numerical analysis of SSME injector atomization and combustion process. 

• Application to bipropellant spray combustion. 

• Numerical simulation of combusting flows with supersonic droplet injection. 
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Table 1. Gas-phase B.C. and droplet-phase I.C. 


Air Mass Row Rate 

355 kg/hr 

Air/Fuel Ratio 

20.17 

Inlet Air Temperature 

310 K 

Droplet Distribution 

Rosln-Rammler 

Sauter Mean Diameter 

127 n m 

Droplet Size Range 

10-290 

Number of Size Range 

15 

Axial Droplet Velocity 

11 m/s 

Tangential Droplet 
Velocity 

6.1 m/s 

Radial Droplet Velocity 

0.5 - 2.5 m/s 

Droplet Temperature 

310 K 
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Table 2. Test Conditions for the Measurement of Hiroyasu and Kadota 
Nozzle diameter: 300 urn Injection Pressure: 9.9 MPa 


Case 

Pgas 

(MPa) 

Pgas 

(kg/m3) 

Vinj 

(m/s) 

Minj 

(kg/s) 

SMD 

(pm) 

1 

1.1 

12.36 

115.80 

0.00688 

42.4 

2 

3.0 

33.70 

102.54 

0.00609 

49.0 

3 

5.0 

56.17 

86.41 

0.00513 

58.8 


Table 3. Test Conditions for the Measurements of Shearer and Groff. 


Pgas 

Pgas 

Vinj 

VOLinj 

Cone Angle 

(kPa) 

( kg/m3) 

(m/s) 

(ml/inj) 

(deg) 

550 

6.36 

60 

0.0165 

60 


Table 4. 

Test Conditions for the Measurement of Yokota et. al. 

Case 

Pinj 

(MPa) 

Pgas 

(MPa) 

Tamb 

(K) 

Minj 

(kg/s) 

Atmosphere 

Evaporating 






Spray 

30 

3.0 

900 

0.00326 

N2 

Burning 






Spray 

30 

3.0 

900 

0.00326 

Air 
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Figure 4.1 Particle dispersion of a nearly-homogeneous flow for 
SSF model(5000 particles) and DDWT model. 
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Figure 4.5 Geometry of the hollow-cone spray combustor(Khalil et. al.) 
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Figure 4.0 Droplet trajectories in kerosene spray Maine fields 
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Figure 4.7 Velocity vectors in kerosene spray flame fields 
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Figure 4.11 Radial profiles of mean axial velocity(S=0.72) 
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Figure 4.12 Radial profiles of mean axial velocity(S=1.98) 
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Figure 4.14 Radial profiles of mean tangential velocity(S=1.98) 
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Figure 4.16 Droplet trajectories in kerosene spray flame fields 



Figure 4.17 Velocity vectors in kerosene 6pray flame fields 
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Figure 4.19 Radial profiles of mean teinperatiire(S=1.2) 




Figure 4.20 Radial profiles of mass fraction of C02(S=1.2) 




Figure 4.21 Radial profiles of mass fraction of 02(S — 1.2) 
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Figure 4.22 Spray parcel distribution in a solid-cone spray(< — 3.0ms) 
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Figure 4.24 Sauter mean diameter versus distance from the injector 
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Figure 4.30 Contours of temperature and oxygen mass fraction in a burning spray 










APPENDIX A. Two-Phase Interaction Source Terms 


The two-phase 
pressed as : 


interaction source terms in the governing equations can be ex- 

NP 

S m , t = £ N p m t ,JdV (A.1) 

P = 1 


iV j 

Sux,i - Y^[N P Tn eViP (v t ) p - iir pdr P *N r (-~)]/dV (A. 2) 

p=l 



Sh,i — L ■+■ ) 



p= i 



- *>,)»/<»• 

(A.4) 

where 

( *i) _«r* , +«i-? +1 +flj 

dt p r 

(A4) 


Here, dV denotes the volume of the computational cell and h p and L are the droplet 
enthalpy and the latent heat of the droplet, respectively. 


To improve the covergence and the numerical stability, the momentum inter- 
action source term, S u ij can be treated implicitly. 


5*+/ = -S u u* +1 + R u (A.5) 

Here, 5„ and R v are obtained by substituting (A. 4) into (A.2): 

. NP 

5, = ^^ N p m p /(At + t p ) (A.6) 

p 

1 NP 

Ru = jy^Npmp/iAt + TpKvf-u'i + F^Tp) (A. 7) 

P 

and m p = \^r z ? p d is the particle mass. The parameters 5„ and R u are momentum 
control volume quantities depending on available particle information at previous 
timestep. 
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APPENDIX B. Numerics of Beta Probability Density Function 

The density- weighted mean mixture properties (<f> m ) at any location evaluated by 
convoluting the property functions with a probability density function, P(£,Xi)’- 






f <f>m(0 P (Z> x i) d t 
J 0 

(Bl) 

where 


C-‘(1 -0‘-‘ 

(52) 


/o'e-'o -?)*-'# 



The denomenator in Eq.(B2) is the Beta function, 5(< z, &). Note that 5(<z, 6) can be 
calculated from the Gamma function, T with the aid of the following relationship. 


B(a, b ) 


r(q)r(6) 
r(a + b ) 


(B 3) 


Substituting Eq.(B2) into Eq.(Bl) yields 

The numerator of Eq.(B5) can be integrated by a trapzoidal rule or Gaussian 
quadrature. However, the significant errors can be produced unless the increments 
in £ are chosen to be quite small. These numerical errors axe due to the spikey 
nature of the Beta pdf when the variance ( g ) of the mixture fraction is small and 
also when the mean mixture £raction(/) is close to 0 or 1. To avoid this problem, 
$ m ( 0 is expressed as polynomials of the conserved scalar(^) . To keep the function 
monotonic, the integration domain (0<£<l)is split into two sections: 0 to 
and ( 3 to 1 , where f , is the stoichiometric value. The property <j> m is expressed as 


n 

when 0 < f < £ s 

(Bo) 

dmn£, 

when £ s < £ < 1 

(B 6) 


where 

77 i = l,...,M (M = index for the property ) 
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n = 1 , N (N — degree of the polynomial ) 


Substituting Eqs.(B5) and (B6) into Eq.(B4) results in 

U. n 

where a and b varying with x, are evaluated from the solutions of the transport 
equations for / and g. The numerator of Eq.(B7) is simplified as Mows: 

T = . - d mn ) [ ‘ r + "‘(l - 


+y,d m n f'r + ‘- l (i-o b -'d( (bs) 

n 

For using IMSL routines, (B8) can be transformed as the convenient expression: 

, j ['r + °-'a -o'-' d( 

T = El*** + (Cmn dm n) 7o * 




= £[d m . + (cmn - d m n)BETAI(s,a + n, b)]BETA(a + n,b) 


(B 9) 


In present study, the degree of the polynomial is used as N = 6. Finally, the mean 
mixture property, ^ m (xj) can be calculated as 




BETA(a,b ) 


(B 10 ) 
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APPENDIX C. Stoichiometric Relations For Hydrocarbon Fuels 

For the hydrocarbon-air mixtures, the irreversible single-step reaction is expressed 
as follows: 

C X H V + (x + |)(0 2 + nJV 2 ) -» xC0 2 + |# 2 O + (i + V - )nN 2 (Cl) 

Here, n is 3.76. In the given reaction process, five species (fuel, 0 2 , N 2 , C0 2 , and 
H 2 0) are participating the mixture composition. Once the mass fraction of fuel 
and oxidizer have been determined from the solutions of the transport equations, 
the mass fraction of the remaining species can be obtained from the following sto- 

ichimetric relations. 



Y"h?o = K 2 (l — K\Yo 2 ~ Yfu) 

(C 2) 


Yco 7 = Kz Yh 7 o 

(C3) 


Y}j . a = 1 - ( Yh 7 o + Yco 7 + Yo 7 + Yfn) 

(C4) 


where 

K '= 1 + n W 0 , 

iW Hi0 

Kl ~ if W„,o + {(x + + zWcoi)} 



xWco 7 
3 " f W H7 o 
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APPENDIX D. Droplet Distribution Models 


The present dilute spray model assumes that the fuel is injected into the combus- 
tion chamber as a fully atomized spray which consists of spherical droplets. The 
droplet-size distribution with the spray is represented by a finite number of size 
ranges. A Nukiyama-Tanasawa distribution and a Rosin-Rammler distribution are 
implemented in the computer code. These distributions have the following forms: 

Nukiyama-Tanasawa Distribution 


dN 

N 


= A( 


SMD 


) 


, — B{D/SMD ) 0 


dD 

SMD 


(D. 1 ) 


where dN and N are the number of computational parcels in the size range from D 
to D+dD and the total number of computational parcels, respectively; SM D is the 
Sauter mean diameter; and a, /3, A, and B are experimental! determined constants. 


Rosin-Rammler Distribution 


dQ gD q ~ l _—( p/x) q 
dD X 9 


(D. 2) 


SMD q ) 

where Q is the fraction of the total volume contained in drops of diameter less than 
D, and X and q are constants. 
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introduction 

A particulate two-phase flow CFD model was developed based 
on the FDNS code (Refs. 1,2,3) which is a pressure based 
predictor plus multi-corrector Navier-Stokes flow solver. 
Turbulence models with compressibility correction (Ref. 4) and the 
wall function models (Ref. 5) were employed as submodels. A 
finite-rate chemistry model (Refs. 6,7) was used for reacting 
flow simulation. For particulate two-phase flow simulations, a 
Eulerian-Lagrangian solution method using an efficient implicit 
particle trajectory integration scheme was developed in this 
study. Effects of particle-gas reaction and particle size change 
to agglomeration or fragmentation were not considered in this 
investigation . 

At the onset of the present study, a two-dimensional version 
of FDNS which had been modified to treat Lagrangian tracking of 
particles (FDNS-2DEL) had already been written and was 
operational. The FDNS-2DEL code was too slow for practical use, 
mainly because it had not been written in a form amenable to 
vectorization on the Cray, nor was the full three-dimensional 
form of FDNS utilized. The specific objective of this study was 
to reorder the calculations into long single arrays for automatic 
vectorization on the Cray and to implement the full three- 
dimensional version of FDNS to produce the FDNS-3DEL code. Since 
the FDNS-2DEL code was slow, a very limited number of test cases 
had been run with it. This study was also intended to increase 
to n umb er of cases simulated to verify and improve, as necessary, 
the particle tracking methodology coded in FDNS. 

Governing Equation 


The gas-phase governing equations of the FDNS module are the 
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Reynolds -averaged Navier-Stokes equations with the addition of 
particle drag forces and heat fluxes in the momentum equations 
and the energy equation, respectively. Due to the effect of 
large density differences between the particles and the 
surrounding gas, the drag force was considered to be the primary 
contribution to the inter-phase momentum exchange. The gas-phase 
governing equations are written as: 

j-^apq/at) = di-pUtf + + s q 

where q - 1, u, v, w, h, k, « and a { for the continuity, 
momentum, energy, turbulence model and chemical species transport 
equations respectively. And, the transformation parameters and 
effective viscosity, /i eff , are given as: 

J * 3(£,»j,r)/a(x,y,z) 
u, - (u/J) (a^j/axj) 
g {J = (a^/ax k ) (a^ j /ax k )/j 
M eff - (/* + /* t )/o q 

The source terms in the governing equations, S q , are given as: 

0 

-P, + ?[M eff ( upj - (2/3)(/, eff 7u) x + Dx 

-Py + V[Mef|(Uj) y ] - (2/3)0i eff VU) y + Dy 

-Pz + *[#»„, (U,),] - (2/3)(^ #ff Vu) I + DZ 

Dp/Dt + h ¥ + H p - u p Dx - v p Dy - w p Dz 

p( P r ~ O 

p( £ /k) t(C,+C 3 P r /c) P r - C 2 €] 


where Dx, Dy and Dz represent the drag forces and n takes on 
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values between 1 and N. u p , v p and w p are the particle velocity 
components. H p is the rate of heat transfer per unit volume to 
the gas phase, hy stands for the viscous heat flux of the gas 
phase. P p stands for the turbulence kinetic energy production 
rate and is written as: 

p r - in t /p) [ (auj/ax, + au/axpVz - 2(au k /ax k ) 2 /3] 

An equation of state, p = p/(RT/MJ, is used to close the above 
system of equations. Turbulent Schmidt and Prandtl numbers, a q , 
for the governing equations and other turbulence model constants 
are given taken from Refs. 4, 6 and 7. 

Finite Rate Chemistry Model 

For gas-phase chemical reaction modeling, a general system 
of chemical reactions is written in terms of the stoichiometric 
coefficients (i/ {j and i^*) and the i-th chemical species name (M ( ) 

of the j-th reaction as 
Z M, - E Hj ' 

\ i 


The net rate of change in the molar concentration of species 
i due to reactions j , X ij; is written as: 

x ij “ C K f j n (/ ,a i/ M wi)‘' 1i " K bj n(pa i /M Hi )*' i ^] 

and the species production rate, (in terms of mass fraction) 

is calculated by summing over all reactions. 

w, ** M wi ZXjj 

J 
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where 

My, - molecular weight of species i 
a, « mass fraction of species i 
p - fluid density 
K f j ■ forward rate of reaction j 

- backward rate of reaction j = K f j/K eJ 
K,j = equilibrium constant 

■ exp{E(f, 'v,j ' - 

f, = Gibbs free energy of species i 
K f * A T 8 exp{-E/RT) 

Finally, the species continuity equations are written as: 


p D t a, - (M«ff/0 Va j] = w i 

where o a (assumed to be 0.9) represents the Schmidt number for 
turbulent diffusion. A penalty function is employed to ensure 
the basic element conservation constraints at the end of every 
time marching step. This is a crucial requirement for the 
numerical stability and accuracy of a CFD combustion model. This 
is accomplished by limiting the allowable changes in species 
concentrations, which are the solutions of the species continuity 
equations, for each time step such that the species mass 
fractions are well bounded within physical limits. The resulting 
limited changes are adjusted so that they are proportional to the 
species source terms. A similar chemistry approach and detailed 
turbulence submodels were reported previously (Ref. 8) . 

Particulate-Phase Equations 

A Eulerian-Lagrangian particle tracking method was employed 
in FDNS to provide effects of momentum and energy exchanges 
between the gas phase and the particle phase. The particle 
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trajectories are calculated using an efficient implicit time 
integration method for several groups of particle sizes by which 
the drag forces and heat fluxes are then coupled with the gas 
phase equations. The equations constitute the particle 
trajectory and temperature history are written as: 


DV,/Dt = (U, - V,)/t d 

Dhp/Dt = Cp, (T -w - T p )/t„ - 6 oit T p V(P p <* p ) 


where U, 
v < 



a 


f 


■ Gas Velocity 

■ Particle Velocity 

= Particle Dynamic Relaxation Time 

- 4 p p dp/ (3 C d p c jU, - V, j) 

■ Particle Enthalpy 

= Particle Heat Capacity 
= Particle Temperature 
*s Gas Recovery Temperature 
= Particle Thermal -Equilibrium Time 
= (P p d p )/ [12 Nu M /(Pr d p )l 

■ Stefan-Boltzmann Constant 

- 4 . 76E-13 BTU/FT 2 -S-R 

* Particle Emissivity » 0.20 — 0.31 
= Radiation Interchange Factor 
= Particle Diameter 

■ Particle Density 


C and Nu stand for drag coefficient and Nusselt number for 
d 

heat transfer which are functions of Reynolds number and relative 
Mach number. Typical correlations are given in Refs. 9 and 10. 
Carlson and Hoglund's correlation (Ref. 9) is written as: 
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C d - (24/Re) (1 + 0.15 Re 0 687 ) (1 + O/ 

[1 + M (3.82 + 1.28 e'^^J/Re] 

Nu - (1 + 0.2295 Re 0 * 55 )/ 

[1 + 3.42 M (2 + 0.459 Re 0 * 55 ) /Re] 

where a - 0.427/M 4 * 63 + 3.0/Re 0 * 88 . A more accurate but more 
complicated correlation for the drag coefficient is provided by 
Henderson (Ref. 10). That is, for Mach a 1, 

C d ■ 24 [Re + S {4.33 + exp(-0.247 Re/S) (3.65 - 1.53 T^T) 

/(I + 0.353 TyT) }]* 1 

+ exp(-0. 5*M/Re 1/Z ) [0.1M 2 + 0.2M 8 + (4.5 + 0.38a) 

/(I + a) ] + 0.6 S [1 - exp (-M/Re) ] 

where S = M(7/2) 1,z is the molecular speed ratio, a = 0.03 Re + 
0.48 Re 1/2 . For Mach £ 1.75, 

C d - [0.9 + 0.34/M 2 + 1 . 86 (M/Re) 1/2 {2 + 2/S 2 

+ 1.058 (T/T) v2 /S - 1/S 4 )] / [1 + 1.86 (M/Re) 1/2 ] 

And, for 1 < Mach <1.75, 

C d “ c d mi + ( 4 /3) (M - 1) (C dN=1>75 - C dM1 ) 

which assumes a linear variation between M * 1 and M = 1.75. 

It has been shown that the Henderson drag law gives better 
motor performance predictions compared with test data. The 
applicability and possible improvement of the Nusselt number 
correlation is currently being actively researched (Ref. 11) . 
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Details of the Particle Solution Method 

In the present two-phase flow model, an independent module 
was employed for the calculation of particle drag forces and heat 
flux contributions to the gas flow field. Subroutines for 
locating the particles and integrating their trajectories are 
called for each particle size group. The drag forces and heat 
fluxes are then saved for every grid point. These forces and 
fluxes are then used to evaluate the particle source terms in the 
gas-phase governing equations. In the present FDNS flow solver, 
two forms of the energy equation (i.e. static enthalpy form or 
total enthalpy form) can be selected. It has been found that 
although either form of energy equation usually gives similar 
solutions, the static enthalpy equation provides better 
definition of the liquid rocket plume shear layers, as shown by 
extensive solutions made for the SSME. A determination of which 
form the energy equation best simulates solid (two-phase) rocket 
motor plumes has not yet been made. 

Particle wall-boundary conditions are treated by using a 
specified fraction of the colliding particles which stick to the 
wall. Particles which stick result in a decreased particle 
velocity normal to the wall for that particle size fraction. 
Therefore, for the particle size fraction which locally collides 
with the wall, part of the particles stick and the other part is 
turned more parallel to the wall. Energy exchange is assumed to 
be due only to the particles which stick. This model of particle 
wall interaction can be improved, but new experimental test data 
must become available in order to do so. 

In the 2-D version of the FDNS flow solver, a fourth-order 
Runge-Kutta method was employed to integrate the particle 
trajectories. After a thorough test of the integration routine, 
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it was found that the explicit scheme can sometimes give diverged 
particle solutions when the source terms become large. 

Therefore, an implicit integration scheme was employed in the 
present model. For convenience, consider the X-component of the 
particle equation of motion. That is, 

dXp/dt - U p 
dUp/dt - A (U c - Up) 

where A - l/t d 

U e - gas velocity 
Up = particle velocity 
X p = particle location 

In finite difference form the above equations can be written as: 

X p <n) - (At/2) [U p <rn1) + Up (n) ] 

U p <n) « AtA [U c - Up (,Hl) ] 

X p <n> + At/2 [U p {n+1> + U p (n> ] 

[U p <n> + AtA U e ]/(l+AtA) 

These two equations are unconditionally stable despite the 
magnitude of the source terms. To provide better time 
resolution, a variable time step size is chosen so that a 
particle would take at least 4 time steps to go across a grid 
cell. 


Y (n*D — 

tt (r*D _ 
P 

or 

Y (n+1) = 
A P 

tj (n*1) m 
P 


The recognition that an improved integration scheme was 
needed for calculating the particle trajectories was a major 
hurdle in developing FDNS-3DEL. The explicit scheme appeared to 
give acceptable solutions, but detailed comparisons to previous 
FDNS-2DEL analyses showed that unacceptable pressure losses were 
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predicted. Several other factors were initially suspected of 
causing this solution behavior. Namely, the turbulence model, 
the form of the energy equation, and the particle drag law were 
initially suspected, and lengthy calculations were made before 
these effects were found not to be the cause of poor results. 
Since the FDNS-2DEL results were found to give good pressure 
field comparisons to conventional nozzle and plume flowfield 
codes (RAMP, SPP, and SPF-II) , the Runge-Kutta method was not 
expected to perform poorly in the FDNS-EL code. Resolving this 
problem consumed much of the resources which otherwise would have 
been used to run a wider variety of test cases. 
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Test cases 

The major test case which was studied was the Tomahawk solid 
rocket motor nozzle analysis, consideration of a plume flowfield 
end of an oxygen-hydrogen coaxial injector was also made. These 
cases are described in the following paragraphs. 

• The Tomahawk Nozzle Flowfield 

The Tomahawk nozzle flowfield was calculated with FDNS-3DEL 
and is shown in Figs. 1-4. This test case was chosen because 
comparable predictions with the FDNS-2DEL and RAMP codes had 
already been performed, and these other solutions were available 
for comparison (Ref. 12). Figures 1-4 show the velocity, Mach 
number, temperature, and water concentration profiles, 
respectively, for the chamber, nozzle, and near plume. The 
chamber flow was approximated to be uniform so that direct 
comparisons with the previous solutions could be made. The FDNS- 
2 DEL solution predicted somewhat lower exit plane centerline gas 
temperatures (2250 °K) than the RAMP solution (2400 °K) . The 
FDNS-3DEL (2470 °K) and RAMP solutions show essentially the same 
exit plane centerline gas temperatures. The raggedness in the 
temperature profile near the centerline in the nozzle appears to 
be due to a weak oblique shock. An apparent non-zero temperature 
normal gradient at the centerline in the subsonic portion of the 
nozzle flowfield is indicated. This is due to a very strong 
effect of the inlet particle flowfield boundary condition. In a 
complete SRM simulation which includes the burning grain, more 
particles would flow down the centerline from the chamber (as 
compared to the uniform flow case) and this subsonic temperature 
contour would probably change shape. The sharp breaks m the 
velocity, Mach number, and temperature contours locate the 
approximate limiting streamline of the particle laded flow with 
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respect to gas only flow which fills the nozzle. Both the static 
and total enthalpy forms of the energy equation give the same 
nozzle solutions. Letting the particles which hit the wall stick 
or elastically reflect give well behaved solutions. The only 
place where there is significant particle impact on the wall is 
at the start of the converging section. The analysis allows 
particles which hit the plume centerline to spectrally reflect, 
in order to account for particles crossing the plume centerline. 
However, particle drag moves the particles very parallel to the 
gas streamlines in the transonic region of the nozzle, such that 
such reflection does not occur in the case being considered. 

• The Tomahawk Plume Flowfield 

The near plume appears to be well predicted with FDNS-3DEL. 
The predicted free shear layer is sharply defined and indicates 
water production from afterburning reactions. Both the static 
and total form of the energy equation were considered. The total 
form of the energy equation indicates a temperature spike at the 
inception of the free shear layer. Better definition of the 
induced flow on the outside of the nozzle would probably 
eliminate such a spike. The static form of the energy equation 
does not exhibit this effect. A Mach number correction to the 
k-« turbulence model was used for this simulation. 

When the Tomahawk plume is calculated for a long distance 
down stream of the exit plane with FDNS-2DEL, excessively rapid 
plume/atmosphere mixing is predicted. This was believed to be 
due to the effect of crossing the Mach disc in the plume and 
thereby creating too much turbulent kinetic energy with the k-« 
turbulence model being used. A similar problem exists when using 
the SPF/II standard JANNAF plume code (Ref. 13). The remedy in 
the SPF/II code is to switch turbulence models between the near 
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and far plumes. An insufficient number of test cases have been 
run with FDNS-3DEL to determine if the Mach number modified 
turbulence model will indeed fix this problem, although the 
solution is better behaved than when the extended k-e turbulence 
model is used. This potential plume prediction problem for far- 
field analysis must be left for future resolution. The FDNS-3DEL 
code should not require any change other than turbulence model 
parameters to adjust the rate of plume/atmosphere mixing. It 
should be noted that the computed results with FDNS-2DEL, at 
first glance, look like the afterburning combustion reaction 
rates are too slow. Actually, so much of the cold atmosphere had 
mixed with the plume that the existence of afterburning was not 
apparent . 

• Liquid Injector Flowfields 

The current version of FDNS-3DEL does not treat mass 
transport from the particle phase to the gas (or continuous) 
phase. Also, the particle phase is treated with a lumped model 
such that the particle temperature is constant throughout the 
particle at any instant of time during the flow through the 
computation field. These restrictions should be removed before 
the code is useful for describing spray combustion. However, the 
spread of a droplet cloud of supercritical fuel or LOX could be 
described with the code without modification, if one is content 
with not describing local mixture ratio changes, i.e. one assumes 
that the supercritical lump remains a lump (or particle) in the 
region of the flow being analyzed. The energy transfer for 
supercritical injection could be easily treated in this manner 
because the heat of vaporization does not have to be considered. 
In fact, models which are based on arbitrarily supplying such 
heat of vaporization (Ref. 14), do not realistically describe 
supercritical spray phenomena. The only reason that such models 
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work at all is that the heat of vaporization evaluated at the 
temperature of the oxygen lump crudely approximates the high heat 
capacity of the liquid-like lump at supercritical pressures. An 
oxygen spray eminating from a single coaxial injector could be 
described with FDNS-3DEL by assuming the oxygen lump to be of a 
constant size and density. A demonstration calculation of this 
nature was considered, but the lump density would be such a very 
strong function of the mean lump temperature that the calculation 
was not performed. If accurate real-gas equation of state models 
were used, the stated oxygen spray simulation would be 
meaningful. Currently, SECA is developing the more general 
property evaluation for a hydrogen-oxygen engine heat transfer 
analysis (Ref. 15) . However, the currently feasible constant 
property analysis was not made, because a reliable two-phase, 3- 
dimensional FDNS-3DEL code was not completed early enough in this 
study . 

Closure 

The calculation of two-phase reacting flows at best is a 
slow process. Several strategies were tried to make this process 
more efficient. Initially, ideal gas flow was computed, then the 
reactions were turned on, and finally the particle trajectories 
were calculated. The entire flowfield was calculated for each of 
these flow conditions. Recently, all of these conditions have 
been treated simultaneously from the beginning of the analysis. 
This procedure works well and results in an overall reduction in 
computation time. For analyzing rocket motors and their 
attendant plumes, it is recommended that the flowfield should be 
broken into subregions for analysis, in order to use the optimum 
step size for the Mach number range within the region. Such a 
restart option has been incorporated in FDNS-3DEL. For example, 
the motor and nozzle should be analyzed first. The computed 


17 



SECA-TR-92-06 


nozzle exit conditions should be used to calculate the near 
plume. The far plume should then be computed. The break between 
the near and far plume should be chosen somewhere between the 
establishment of the complex near field shock structure and the 
essentially balanced jet, predominately mixing dominated far 
field. The development of a parabolized version FDNS-3DEL to 
initially predict large plume structures and other large 
flowfields is also recommended. 
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Conclusions 

1. A two-phase, finite-rate CFD code (FDNS-3DEL) was developed 
and vectorized. The Tomahawk nozzle test case indicates the 
CFD solution accurately simulates this flow. 

2. Particle mass transfer effects are not currently included in 
the current code. The inclusion of these effects would be 
relatively simple. 

3. More test cases should be run to establish the range of 
validity of the calculation procedure. The mechanics of the 
Euler-Lagrange calculation appear to be in good working 
order. Secondary effects, such as turbulent-mixing/ shock- 
structure interaction require further study with more test 
cases. However, it should be noted that suitable 
experimental data to verify many of these complex flow 
interactions are not now available. The best one can 
currently do is compare CFD solutions to SPF-II type 
analyses. 

4. Analyzing large, complex flowfields with any two-phase, 

finite-rate CFD code is a time consuming process, therefore 
utilization of all methods which would expedite such 
analyses should be considered. Analyzing the flowfields 
with carefully selected subregions and developing 
parabolized versions of the CFD codes are two such 
computational aids which should be employed. 
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